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>^: 

' Many astrophysical phenomena are highly subsonic, requiring speciahzed nu- 

merical methods suitable for long-time integration. In a series of earlier papers we 
described the development of MAESTRO, a low Mach number stellar hydrodynam- 
ics code that can be used to simulate long-time, low-speed flows that would be 
prohibitively expensive to model using traditional compressible codes. MAESTRO 
(-H ' is based on an equation set derived using low Mach number asymptotics; this 

Qh' equation set does not explicitly track acoustic waves and thus allows a significant 

O ■ increase in the time step. MAESTRO is suitable for two- and three-dimensional 

^ • local atmospheric flows as well as three-dimensional full-star flows. Here, we 

^ ■ continue the development of MAESTRO by incorporating adaptive mesh refinement 

(AMR). The primary difference between MAESTRO and other structured grid AMR 
^ ' approaches for incompressible and low Mach number flows is the presence of the 

time-dependent base state, whose evolution is coupled to the evolution of the 
full solution. We also describe how to incorporate the expansion of the base 
state for full-star flows, which involves a novel mapping technique between the 
Q ■ one-dimensional base state and the Cartesian grid, as well as a number of overall 

P ■ improvements to the algorithm. We examine the efficiency and accuracy of our 

adaptive code, and demonstrate that it is suitable for further study of our initial 
scientific application, the convective phase of Type la supernovae. 
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Introduction 



Many astrophysical phenomena of interest occur in the low Mach number regime, where 
the characteristic fluid velocity is small compared to the speed of so und. Some well-known 



examples are the convective phase of Type la superno yae (SN la) (iHofiich fc SteinI 12002 



Kuhlen et al. 20061: Zingale et al.ll2009l ). classical novae (IGlasner et al.l 



20071) 



.^ , , V— i ^^^^^^, , convection in 

stars ( Meakin &: ArnettI 120071 ). and Type I X-ray bursts ( Lin et al.l 120061 ). Such problems 



require a numerical approach capable of resolving phenomena over time scales much longer 
than the characteristic time required for an acoustic wave to prop agate across the com- 
putational domain. In a series of papers fsee lAlmgren et al.l (l2006al) — henceforth Paper I, 
Almgren et al. ( 2006bl)— henceforth Paper II, lAlmgren et al.l ( 20081 ) — henceforth Paper III, 



and IZingale et al.l ( 120091 ) — henceforth Paper IV), we have described the initial development 
of MAESTRO, a low Mach number hydrodynamics code for computing stellar flows using a 
time step constraint based on the fluid velocity rather than the sound speed. MAESTRO is 
suitable for two- and three-dimensional local atmospheric flows as well as three-dimensional 
full-star flows. All simulations are performed in a Cartesian grid framework, but rely on the 
presence of a one-dimensional radial base state that describes the average state of the star or 
atmosphere. Starting with the development of the low Mach number equation set (see Paper 
I), we demonstrated how to capture the expansion of the base state in a local atmospheric 
simulation in response to large-scale heating (Paper II) and incorporate reaction networks 
(Paper III). In Paper IV, we presented the initial application of MAESTRO, following the last 
two hours of convection inside a white dwarf leading up to the ignition of a SN la using a 
three-dimensional, full-star simulation with a base state that is constant in time. 

In general, astrophysical flows are highly turbulent. In the ca se of the convective p eriod 
preceding a SN la explosion, the Reynolds number is (9(10^^) (jWoosley et al.ll2004l ). far 
larger than can be modeled on today's supercomputers. Nevertheless, to understand the 
role of turbulence in these events, we must use increasingly more accurate simulations. In 
this paper we describe how to incorporate adaptive mesh refinement (AMR), in which we 
locally refine the Cartesian grid in regions of interest, to allow us to efficiently push to higher 
spatial resolutions and better capture the turbulent flow in critical regions of the simulation. 
The primary difference between MAESTRO and other structured grid AMR approaches for 
incompressible and low Mach number flows is the presence of the time-dependent base state, 
whose evolution is coupled to the evolution of the full solution. We also describe how 
to incorporate a time-dependent base state for full-star problems, which involves a novel 
mapping technique between the one-dimensional base state and the Cartesian grid. This 
allows us to properly capture the effects of an expanding base state in full-star simulations. 
We have also made a number of overall improvements to the algorithm, and all together, 
these enhancements will allow us to compute more efficient and accurate solutions for our 
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target applications, including the convective phase of SNe la and Type I X-ray bursts. 

This paper is divided into several sections, along with three detailed appendices. In §21 
we present the governing equations. In §31 we give an overview of the methodology, referring 
the reader to the appendices for full details. In §11 we describe the new mapping procedure 
between one-dimensional and three-dimensional data structures in full-star problems. In §5l 
we discuss the extension of the algorithm to include AMR. In §6l we describe the results of our 
test problems. We conclude with §3 which includes future plans for scientific investigation. 



2. Governing Equations 

Stellar flows are well characterized by the compressible Euler equations (i.e., viscosity 
effects are negligible). These equations model all compressibility effects in a fluid, and allow 
for the formation and propagation of shocks. For low speed convective flows in a hydrostat- 
ically stratified star or atmosphere, we do not need to explicitly follow the propagation of 
sound waves. However, we do need to include large-scale compressibility effects such as the 
expansion/contraction of a fluid element as it changes altitude in the stratified background, 
and the local changes to the density of the fluid element through heating and compositional 
changes. By reformulating the equations of hydrodynamics to filter out sound waves but 
preserve the correct large-scale fluid motions and hydrostatic balance, we can retain the 
compressibility effects we desire while allowing for much larger time steps than a corre- 
sponding compressible code. The full derivation of the low Mach number hydrodynamics 
equations is given in papers I through III. The resulting equations are: 

^Mtl = -V ■ (pxtv) + pLu,, (1) 

^ . _U.VU-iv.-^^^se„ (2) 

ot P P 

^ = -V-(p/^U) + ^ + pi7nuc + pi^e.t, (3) 

where p, U, and h are the mass density, velocity and specific enthalpy, respectively, and 
Xk are the mass fractions of species k with associated production rate Uk- The species are 
constrained such that Ylk-^k = 1 giving p = X]fc(P^fc) ^^^ 

| = -V.(PU). (4) 

The source terms Hcxt and H^uc are the external heating rate and nuclear energy generation 
rate per unit mass. The pressure is decomposed into a hydrostatic base state pressure. 
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Po = Po{r,t), and a dynamic pressure, n = 7r(x, t), such that |7r|/j)o = 0{M'^) (we use x to 
represent the Cartesian coordinate directions of the full state and r to represent the radial 
coordinate direction for the base state). We also define a base state density, po = Poij'^'t), 
which is in hydrostatic equilibrium with pq, i.e., Vpo = —Po9^r, where g = g{r,t) is the 
magnitude of the gravitational acceleration and e^ is the unit vector in the outward radial 
direction. 

Mathematically, this system must still be closed by the equation of state which we 
express as a divergence constraint on the velocity field (see Paper III), 

V.(AU)^/,.(s-=i-f). (5) 

where Po is a density-like variable that carries background stratification, defined as 

/3o(r, t) = po(0, t) exp (^£ ff;^^^') ' (6) 

and Fi the lateral average (see §4.ip of Fi = (i(logp)/(i(logp) at constant entropy. The ex- 
pansion term, S, incorporates local compressibility effects due to heat release from reactions, 
compositional changes, and external sources, 

S = -ay~] ^kUJk H y~] PXk^k + o-i^rmc + o-i^cxt- (7) 

k PPfk 

where pxfe = dp/dXk\p^T,x,,,^,^ ^k = dh/dXk\p^^x^^^^ ,Pp = dp/dp^^^^, and a = Pr/ipCpPp), 
with pt = dp/dT\ -^^ and Cp = dh/dT\ -^^ is the specific heat at constant pressure. 

It is important to note that if the Mach number of the fluid in a numerical simulation 
becomes 0(1), through large acceleration due to buoyancy or nuclear energy generation, for 
example, the solution of these equations would no longer be physically meaningful. The low 
Mach number equations do not enforce that the Mach number remain small; rather, if the 
dynamics of the flow are such the Mach number does remain small, then these equations are 
valid approximations for the evolution of the flow. 

As in Papers II and III, we decompose the full velocity field into a base state velocity, 
Wo, that governs the base state dynamics, and a local velocity, U, that governs the local 
dynamics, i.e., 

U = wo(r,t)e, + U(x,t). (8) 



with (U ■ e^) = and wq = i\] • e^). The velocity evolution equations are then 

dwQ dwo 1 diTo 

-^- = -^0^ ^-, (9) 

ot or po Or 

— - = _u ■ VU - U • e J — — e^ Vtt H --e^ ger. (10) 

at \ J or p pq or p 
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where ttq is the base state component of the perturbational pressure. By laterally averaging 
to equation (|5]), we obtain a divergence constraint for wq: 

V ■ i^oWoBr) = /3o (S - =i-^) . (11) 

The divergence constraint for U can be found by subtracting (ITT|) into (|5]), resulting in 

W-((3o\j)=(3o{S-S). (12) 



In the present paper, we revert back to the method introduced in paper II and define 
a base state enthalpy, (p/i)o- We use po and {ph)o to define the perturbational quantities 
p' = p — Po and (ph)' = (ph) — {ph)^, which are predicted to the Cartesian edges to compute 
fluxes for the conservative updates of po and {ph)o. Experience has shown that by advancing 
perturbational quantities, the slope limiters are more effective at reducing numerical oscil- 
lations since they are being applied to a normalized signal, rather than a signal that spans 
many orders of magnitude over a small number of cells. This is a departure from paper III 
where we predicted temperature to the Cartesian edges. Evolution equations for po and 
(p/i)o are designed so that po and (p/i)o will remain the average over a layer of constant 
radius of p and (ph). The fluxes for (pXk) are computed by first predicting po, p', and X^ to 
time-centered Cartesian edges. The flux for (p/i) is computed by first predicting (p/i)o and 
(ph)' to time-centered Cartesian edges. 

We now derive the equations used to predict the time-centered Cartesian edge values in 

the actual algorithm. The species evolution equation is found by combining equations ([T]) 

and dH): 

dX 

-^ = -V-VXk + oOk. (13) 

at 

The base state evolution equations for density and enthalpy can be found by averaging (j4]) 

and ([3]) respectively over a layer of constant radius, resulting in 

^ = -V-(po«;oe.), (14) 

d{ph)o 



dt 



-V ■ [{ph)oWoer] + ^ + pHnuc + P^cxt- (15) 



where ip is the Lagrangian change in the base state pressure defined a.s ip = Dopo/Dt = 
dpo/dt + wodpo/dr and is related to the total pressure by 

^ = ^ + U-Vpo. (16) 
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Subtracting the base state evolution equations from the corresponding full state equations 
yields 



'dt 

djph)' 
dt 



-U ■ Vp' - p'V ■ U - V • (poUj , (17) 

-U ■ V(p/i)' - (p/i)'V ■ U - V • \{ph)oij] + U ■ Vpo 



+ {pHnnc - pHnuc) + (P-f^cxt - pHc 



In our treatment of enthalpy, we split the reactions and external heating from the hydrody- 
namics, i.e., during the hydrodynamics step, we neglect the Uk, -f^nuc, and H^xt terms. Also, 
in our treatment of species, we similarly split the reactions from the hydrodynamics. 

While equation (TTll) properly captures the change in po due to atmospheric expansion 
caused by heating, it neglects changes that can occur due to significant convective overturn- 
ing. We impose the constraint that p' = for all time. In Paper III, we quantified the drift 
in p' by introducing rjp in the equation 



dfy_ 
'dt 



-V ■ {r],er). (19) 



However, we incorrectly derived r]p by assuming V ■ (p'^o^r) = 0, when in general this is 
not true since p', when predicted to time-centered edges, does not in general satisfy p' = 0. 
Therefore, the correct expression is r]p = (p'TJ-e^). In practice, we correct the drift by 
simply setting po = p after the advective update of p. However we still need to explicitly 
compute rip since it appears in other equations. 



3. Overview of Numerical Methodology 

We shall refer to local atmospheric flows in two and three dimensions as problems 
in "planar" geometry, and full-star flows in three dimensions as problems in "spherical" 
geometry. The solution in both cases consists of the Cartesian grid solution (U, p, h,Xk,T) 
and the one- dimensional base state solution {wq, po, {ph)o,po), all of which are cell-centered 
except for wq, which is edge-centered. Edge-centered data is denoted by a half- integer 
subscript. See Figure [T]for a representation of each grid structure. The time step index is 
denoted as a superscript. 

For planar problems, e,, is in alignment with the Cartesian grid unit vector in the 
outward radial direction, e^ (in two dimensions) or e^ (in three dimensions). We choose 
Ar = Ax so that there will be a simple, direct mapping between the radial array and the 
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Fig. 1. — (Left) For data on the Cartesian grid (shown here in two dimensions), we use a 
cell-centered convention with indices i,j, k (in three dimensions). Edges are denoted with a 
half-integer. (Right) The base state lives on a radial array and uses a cell-centered convention 
with index j. Edges are denoted with a half-integer. 

Cartesian grid. For spherical problems, e^ is not in alignment with any Cartesian coordinate 
direction. Our choice of Ar can be independent of Ax; as in Paper IV, we use 5Ar = 
Ax. Note that for spherical problems, we place the center of the star at the center of the 
computational domain, and therefore the center of the star lies at a corner where 8 Cartesian 
cells meet. See Figure [2] for an illustration of the relationships between the radial array and 
the Cartesian grid for spherical and planar geometries. 

The time- advancement algorithm uses a predictor-corrector formalism. In the predictor 
step, we compute an estimate of the expansion of the base state, then compute a preliminary 
estimate of the state at the new time level. In the corrector step, we use this preliminary 
state to compute a new estimate of the expansion of the base state, and then compute 
the final state at the new time level. We incorporate reactions and external heating using 
Strang-splitting. As in previous papers, our algorithm is second-order in space and time. 

The full details of the algorithm are presented in Appendix [Al The main algorithm 
description in Appendix IA.4I is similar to the description in Paper III, but has been signifi- 
cantly updated to show how we incorporate the time-dependent spherical base state. There 
are numerous other improvements we have made to the algorithm since papers III and IV, 
which are described in Appendix lA.li Note that these changes have also been incorporated 
into the main algorithm description. Overall: 



Appendix lA.ll is a summary of algorithmic changes since papers III and IV. 

Appendix IA.2I describes how we compute and discretize gravity. 

Appendix IA.3l is a description of shorthand notation we use in describing the algorithm. 




Fig. 2. — (Left) For problems in spherical geometry, there is no direct alignment between the 
radial array cell centers and the Cartesian grid cell centers. (Right) For problems in planar 
geometry, there is a direct alignment between the radial array cell centers and the Cartesian 
grid cell centers. 

• Appendix IA.4I steps through the algorithm in detail. 

• App endix I A . 5 1 descr ib es special treatment given to low density regions in the simulation. 



4. Mapping 

At many points in the algorithm, we need to map the full state on the Cartesian grid onto 
a one-dimensional radial array, and vice-versa. Since Paper IV, we have greatly increased the 
accuracy of the numerical mapping to and from these data structures for spherical problems, 
most notably the lateral average routine described below. We refer to the procedure for 
mapping a cell-centered Cartesian field to a cell-centered radial array as a "lateral average" , 
and we refer to the procedures for mapping an edge- or cell-centered radial array to an edge- 
or cell-centered Cartesian grid as a "fill" . 



4.1. Lateral Average 

For any Cartesian cell-centered field, 0, we define =Avg(0) as the lateral average 
over a layer at constant radius r, Qh, as 






dn. (20) 
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planar: This is a straightforward arithmetic average of cells at a particular height since the 
radial cell centers are in alignment with the Cartesian grid cell centers. 

spherical: It can be shown that any Cartesian cell center is a radius f^ = Axy3^+2m 
from the center of the star, where m > is an integer. For example, the Cartesian cell 
with coordinates (i, j, k) = (1, 1, 1) relative to the center of the star lies at a distance 
of Ax^/^/4 + 6 from the center of the star, corresponding to ?Ti = 3. The Cartesian 
cells with coordinates {i,j,k) = (2, 0, 0), (0, 2, 0), or (0,0,2) relative to the center of 
the star also lie at that same distance. For the 384^ resolution examples in this paper, 
we have verified that a non-zero set of Cartesian cell centers map into each radius f^ 
until m is large enough to correspond to a radius larger than half the width of the 
computational domain (i.e., the edge of the domain, not the corner of the domain). 
Figure |3] shows the number of Cartesian cells that map into each radius fm, which we 
refer to as the "hit count" , for a 384^ domain. We use this mapping to help construct 
the lateral average, using the following steps: 

1. Create an itemized list, 0^, where each element is associated with a radius f^ = 



Axy'^/i + 2m from the center of the star. 

2. For each 0^, compute the arithmetic average value of the Cartesian cells whose 
centers lie at the associated radius. As an additional element in the itemized list, 
include the center of the star (corresponding to a radius of r = 0). Compute this 
additional value of at this location using quadratic interpolation with (pQ, 0i, and 
a homogeneous Neumann condition at r = as the stencil points. Note that for 
very large values of m, it is possible that no Cartesian cell centers exist at a radius 
fm (i.e., the hit count is zero). If so, we say that (pm has an undefined/invalid 
value, and we ignore such values for the rest of this procedure. 

3. To compute the lateral average, use quadratic interpolation using the value in the 
itemized list with the closest associated radius, (pk, and the nearest values above 
and below, 4)k+ and 0a:_, using divided differences: 

4>k.-4>k 4>k-<f>k_ 



(p{r) = (pk- + — [r - TkJ H 2 [r - rk_){r - r^), 

where fk_ , ffc, and f^ are the three radii associated with 0fc_ , 0^, and 0^ . Finally, 



constrain 0(r) to lie within the range of 0fc_, 0^, and 0^^ so as to not introduce 
any new maxima or minima. 

In §6.11 we show the improvement of this averaging procedure over the Paper IV 
procedure. 



-10- 



6000 



5000 



4000 



3000 



2000 



1000 








10000 



20000 30000 40000 
itemized list index, 'm' 



50000 



60000 



Fig. 3. — The number of Cartesian cells whose centers lie at a radius fm (i.e., the "hit 
count") for a 384^ domain vs. the itemized list index, m. Indices m < 18,432 correspond to 
locations within half the width of the computational domain. A non-zero set of Cartesian 
cell centers maps into the radius associated with every m < 37, 912, which corresponds to 
approximately 0.72 times the width of the computational domain. The inset plot is a zoom-in 
of the innermost 75 values of m. 

4.2. Fill 

There are four different mappings from a one-dimensional radial array to the three- 
dimensional Cartesian grid; below we describe the procedures for planar and spherical ge- 
ometries separately. 

planar: 

1. To map a cell-centered radial array onto Cartesian cell centers, we use direct-injection 
since the radial cell centers are in alignment with the Cartesian cell centers. 

2. To map an edge-centered radial array onto Cartesian cell centers, we average the two 
nearest radial edge-centered values. 



3. To map a cell-centered radial array onto Cartesian edges with normal in the radial 
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direction, we use 4th order spatial interpolation. For example, in two dimensions, 

7 1 

^i,i+y2 = Y^(<^i + </'i+i) - Y^(<^i-i + <^i+2)- (21) 

We constrain (pij+y^ to lie between the interpolated values, and lower the order of 
interpolation near domain boundaries. For the Cartesian edges transverse to the base 
state direction, we use direct-injection since the radial cell centers are in alignment 
with these Cartesian edges. 

4. To map an edge-centered radial array onto Cartesian edges, we use direct-injection 
on Cartesian edges normal to the base state direction since the radial edges are in 
alignment with these Cartesian edges. For the remaining Cartesian edges, we average 
the two nearest radial edge-centered values. 

spherical: 

1. To map a cell-centered radial array onto Cartesian cell centers, we use quadratic inter- 
polation from the nearest three radial cell centers (see Figure H^). This is a departure 
from Paper IV, in which we used piecewise constant interpolation. 

2. To map an edge-centered radial array onto Cartesian cell centers, we use linear inter- 
polation from the nearest two points (see Figure Hb). 

3. To map a cell-centered radial array onto Cartesian edges, we first map the radial array 
onto Cartesian cell centers (see 1.), then average the two neighboring centers to obtain 
the Cartesian edge values (see Figure Ht). 

4. To map an edge-centered radial array onto Cartesian edges, we first map the radial 
array onto Cartesian cell centers (see 2.), then average the two neighboring centers to 
obtain the Cartesian edge values (see Figure Ht). 



5. Adaptive Mesh Refinement 

Our approach to AMR uses a nested hierarchy of logically rectangular grids with suc- 
cessively f iner grids at high e r leve ls. This is based on the strategy introduced for gas dy- 



na mics bv iBerger fc Co 



bvlAlmgren et al.l (11998 



ellal ( 119891 ) ■ extended to the incompressible Navier-St qkes equations 



) , and extended to low Mach number reacting flows by lPember et al. 



( 119981 ) and iDay fc Belli (J2000[ ). We refer the reader to these works for more details. The 
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Fig. 4. — Illustrations of the fill operation for spherical geometry, (a) To fill the Cartesian 
cell center (fill type 1), represented by the square, from radial cell-centered data, represented 
by the circles, we use quadratic interpolation from the nearest three points, (b) To fill the 
Cartesian cell center (fill type 2), represented by the square, from radial edge-centered data, 
represented by the circles, we use linear interpolation from the nearest two points, (c) To 
fill a Cartesian edge (fill types 3 and 4), represented by the squares, first fill the Cartesian 
cell centers, represented by the circles, then average the two neighboring cell centers. 

key difference between our method and these earlier methods stems from the presence of 
a one-dimensional base state whose time evolution is coupled to that of the full solution. 
To the best of our knowledge, there are no existing AMR algorithms for astrophysics or 
any other field, for flows with a time- dependent base state coupled to the full solution. For 
simplicity, we present a version of the algorithm with no subcycling in time, i.e., the solution 
at all levels is advanced with the same time step. 

We first summarize our AMR approach without the base state, then discuss how the 
base state is incorporated in both the planar and spherical cases. 



5.1. Creating and Managing the Grid Hierarchy 

At each time step the state data is defined on a nested hierarchy of grids, ranging from 
the base level (£ = 1), which covers the entire computational domain, to the finest level 
{i = ^max)- At each level there is a union of non-intersecting rectangular grids with the same 
spatial resolution. For simplicity, we require that the cells composing the grids be square 
[Ax = Ay = Az), and that the refinement ratio between levels be 2. The grids in the 
interior of the computational domain are required to be properly nested, i.e., the union of 
grids at level i+1 is completely contained in the union of grids at level i. Additionally, in the 
interior, we require that each grid at level £+ 1 be a distance of at least two level i cells from 
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the boundary between level £ and level i — 1 grids; this allows us to always fill "ghost cells" 
at level £ + 1 from the level i data (or the physical boundary conditions, if appropria t e). W e 



initialize the grid hierarchy and regrid following the procedure outlined in lBell et al.l ( 1l994l ). 



A user-specified error estimation routine is used to tag cells whe re more resolution is d esired 



The tagged cells are grouped into rectangular patches following iBerger fc RigoutsosI (jl99ll ). 
and subsequently refined to create new grids at next level. Refinement continues until the 
maximum level is reached. 

During Step Op grids at all levels are filled directly from the initial data. As the 
simulation progresses, we periodically check our refinement criteria and regrid as necessary. 
This regridding takes places during Step 12, before computing the next time step. Newly 
created grids are filled by using data from previous grids at the same refinement level (if 
available) or by interpolating from underlying coarser grids. 



5.2. Communication between levels 

Since we use the same time step to advance the solution at all levels, much of the 
com plication associated w ith synchronization of data between levels in a subcycling algorithm 



see 



Almgren et al.lll998l ) is eliminated. The MAC projections in Step 3 and 7 enforce that 
-jjADV,* gj^^ ijADV^ respectively, on any coarse edge underlying fine edges are the average of 
the values on the fine edges. Similarly, the nodal projection in Step 12 enforces that at any 
coarse node underlying a fine node, the value of on the coarse node is identical to the value 
on the fine node above it. The additional communication of data between levels occurs as 
follows: 

• Before any explicit operation at level £ > 1, data in ghost cells at that level are filled by 
interpolating from level i—1, or imposing physical boundary conditions, as appropriate. 

• Edge-based fluxes at level £ < imax that underlie edges at level i + 1 are defined to be 
the average of the fluxes on level i + 1 a.t that edge. This enforces conservation. 

• After any update to the solution, data at finer levels is conservatively averaged onto 
the underlying coarse grid cells, starting at the finest level. 



^The "Step" notation is used in describing the full details of the algorithm in Appendix lA. 41 
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Fig. 5. — (Left) For multi-level problems in planar geometry, we force a direct alignment 
between the radial array cell centers and the Cartesian grid cell centers by allowing the 
radial base state spacing to change with space and time. (Right) For multi-level problems in 
spherical geometry, since there is no direct alignment between the radial array cell centers 
and the Cartesian grid cell centers, we choose to fix the radial base state spacing across 
levels. 



5.3. AMR with a Time-Dependent Base State 

Our specific treatment of AMR is guided by our initial scientific applications, including 
Type I X-ray bursts and the convective phase of SNe la, as well as numerical concerns, most 
notably the presence of the one- dimensional base state. Our treatment of the base state in 
an AMR framework differs for planar vs. spherical problems. 

For planar problems, our approach is to define a radial base state array with variable 
mesh spacing. A general localized fine Cartesian grid would require either a base state that 
exists at multiple resolutions at a particular height, or an interpolation algorithm to obtain 
the base state value at a particular height if Ar ^ Ax. Both of these methods pose problems, 
as they generate oscillations in perturbational quantities (such as p', {pK)' and the S — S term 
on the right hand side of the divergence constraint) since the lateral average routine is only 
defined when the base state is aligned with the Cartesian grid across the width of the domain. 
Any attempt at interpolation will cause oscillations in the perturbational quantities directly 
related to the interpolation error. We have found that such oscillations can be detrimental 
to the results. With these issues in mind, we choose to only allow fine grids to exist that 
span the width of the domain. This way, the base state exists as a single seamless entity 
with multiple resolutions depending on height (see Figure E]). We will take advantage of this 
type of grid structure in our studies of Type I X-ray bursts. 

Next, we define ghost cell values for the finer base state levels, and fill these values 
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by interpolating coarser data. This makes the algorithm directly compatible with the one- 
dimensional time-centered edge state calculation used in Advect Base Density g and 
Advect Base Enthalpy. In particular, the slopes can be used with a consistent stencil at 
each level, that is not dependent on the data from any other level once the ghost cells are 
set. 

Finally, whenever we regrid the Cartesian grid data, we regrid the base state to match 
the grid structure of the Cartesian grid. Then, we set Po = P and compute po using Enforce 
HSE. To compute ip and Wq on the new base state array, we use piecewise linear interpolation 
of the coarser data to fill any new fine radial cells/edges. 

For spherical geometry, we first note that even in the single-level case the radial base 
state is not aligned with the Cartesian grid. Therefore, we use a base state with a fixed Ar 
for all levels (see Figure EJ. As in the single-level algorithm, we choose Ar = Ax/S, but 
here. Ax corresponds to resolution of the Cartesian grid at the finest level. 

Our next consideration is defining the radial average. First, we first create an itemized 
list associated with each level of refinement using only Cartesian cells that are not covered by 
cells at a finer level. At this point one option would be to merge the lists and proceed as in 
the single-level algorithm; this was tested and found to be problematic. Instead, we choose 
the list from a chosen particular level and define the average using quadratic interpolation 
with only this list, as in the single-level case. To decide which list to use, we first examine 
the three points that would be used by quadratic interpolation at each level. The guiding 
principle is to avoid using interpolation points with low hit counts. Thus, at each level we 
find the minimum hit count of the three points; the level which has the largest minimum hit 
count is the level whose list we use for interpolation. We note that this multi-level average 
works particularly well when the center of the star is fully refined. This is the case since 
near the center of the star, there are relatively few Cartesian cells that contribute to each 
radial bin, so by fully refining the center of the star, we ensure that the multi-level averaging 
procedure retains the accuracy of a single-level spherical average near the center. For our 
studies of SNe la, we will take advantage of this fact by always refining the center of the star, 
which is our region of interest (see Figure [12] in §6.6^ as an example). In §6.11 we present 
numerical tests of the new multi-level averaging procedure for spherical geometry. 

For both planar and spherical problems, after regridding the Cartesian grid, we make 
the state thermodynamically consistent by computing T = T{p, h, X^) (for planar problems) 
or T = T{p,pQ,Xk) (for spherical problems). Then, we recompute Fi and /3o as described in 
Steps 10 and 11. 



^Thc boldface notation refers to numerical modules we have described in Appendix lA. 31 
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5.4. Parallel Implementation 

We parallelize the algorithm by distributing the grids on each level across processors. 
Each grid carries a perimeter of ghost cells that are filled from neighboring grids at the same 
level or interpolated from coarser grids as needed. This allows the data on each grid to be 
updated independently of the other grids. A typical grid is large enough (e.g., 32^ cells) 
that for explicit operations the cost of computation within each grid greatly exceeds the cost 
of communication between grids. The linear solves necessary for the MAC projection and 
approximation nodal projection have higher communication costs, but we still obtain good 
parallel efficiency for the overall algorithm. A scaling study for MAESTRO can be found in 



Almgren et all (120071 ). 



Since the one-dimensional base state arrays are so much smaller than the three-dimensional 
arrays holding the full solution, each processor owns a copy of the entire one- dimensional 
base state arrays. Operations such as averaging to define base state quantities require a 
collection operation among grids, followed by a distribution of the average state to each 
processor. 



6. Test Problems 

We have developed a suite of test problems in order to test various aspects of our code. 

» In %.1\ we show that our new mapping procedure from ^is much more accurate than 
the mapping procedure from Paper IV. 

» In §6.2[ we show that we are able to properly capture the expansion of the base state 
in a three-dimensional full star simulation due to heating at the center of the star. 

» In §6.3[ we show that our multi-level algorithm is second-order accurate in space and 
time by tracking a hot bubble rising in a white-dwarf environment. 

» In §6.4[ we show that an adaptive algorithm in three-dimensional planar geometry can 
properly track a hot bubble rising in a white-dwarf environment. 

» In §6.5[ we demonstrate that a multi-level, two-dimensional planar simulation will 
properly capture the expansion of the base state due to a heating layer, and also that 
a multi-level simulation is able to capture the same fine-scale structure as a single-level 
simulation at the same effective resolution over a short time. 
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• In §6.6[ we demonstrate that a full-star simulation with AMR can be used to study 
the dynamics of convection in white dwarfs. 

For test problems §6.2j - §6.6l we compute flows in which the density spans at least four 
orders of magnitude. The large drop in density in the upper atmosphere results in high 
velocities due to conservation of momentum. This region should not affect the dynamics 
below the surface in the convecting regions of the star. However, because the time step in 
the low Mach number code is limited by the highest velocity in the computational domain, 
the efficiency gains of the low Mach number algorithm are reduced if those velocities persist. 
We employ a sponging technique to damp such velocities. Da mping te c hniqu es are commonly 



used in modeling atmospheric convection (see, for example iDurranI fll990l )). In Paper IV, 
for full-star convection, we explored the effects of sponging the velocity beginning at two 
different heights to demonstrate that the dynamics in the upper atmosphere do not affect 
the convecting regions of the star. 

Full details for the sponge implementation can be found in papers III and IV, but in 
summary, we add a forcing term to the velocity update before the final projection. We use 
the parameters rgp, r^d, and k to describe the sponge. The sponge forcing turns on at radius 
Tsp and reaches half of its peak strength at radius r^^- We can control the strength of the 
forcing with the parameter, k. For full-star problems, we also use an outer sponge which 
prevents the velocities near the domain boundaries from becoming too large. 

For all of these tests, we use a public ly available, general stellar equation of state 



(JFryxell et al.ll2000l : iTimmes &: Swestyll2000l ). with contributions from ions, radiation, and 



degenerate/relativistic electrons. 



6.1. Mapping 

To test the new spherical fill and lateral average routines from §H we first create a unit 
cube with 384^ resolution and no refinement. We create a radial array with Ar = Ax/S, 
and initialize the radial array cell centers with the Gaussian profile ^exactl'") = e~^°^' . We 
map 0cxact to Cartesian cell centers with the fill operation, then compute the lateral average 
of this Cartesian field, Avg(0). We repeat this process by choosing the grid with two levels 
of refinement used in the full-star simulation test in §6.6[ shown in Figure | 



Figure El (left) shows the relative error between 0exact and Avg(0) for the new mapping 
procedures and the mapping procedures from Paper IV for the single- level test. The new 
mapping procedure greatly decreases the relative error. Figure [6] (right) is a zoom-in of the 
relative error for the new mapping. For the single-level grid, the relative error is (9(10~®) for 
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Fig. 6. — (Left) Relative error between (/)oxact and Avg(</)) from averaging a mapped Gaussian 
profile centered on a unit cube with 384^ resolution and no refinement for the current mapping 
(solid line) and the mapping from Paper IV (dashed line). The new mapping procedure 
greatly decreases the relative error. (Right) A zoom-in of the relative error for the new 
mapping. The red markers correspond to the test with no refinement. The relative error 
is O(10~^) for r G [0,0.036] and is 0(10"^^) for r G [0.036,0.7]. The green dots shows the 
relative error for a test with two levels of refinement using the grid structure shown in Figure 
[T2l The relative error is (9(10~^) for r G [0,0.7], which is still a vast improvement when 
compared to the Paper IV mapping applied to a single-level simulation. 



r G [0, 0.036] and is (9(10~^^) for r G [0.036, 0.7]. For the test with two levels of refinement, 
the relative error is (9(10~^) for r G [0, 0.7], which is still a vast improvement when compared 
to the Paper IV mapping applied to a single-level simulation. 



6.2. Spherical Base State 

To test the base state expansion for spherical geometry, we perform a series of tests 
similar to those in Paper II which tested the base state expansion in planar problems. We 
run the same test using thre e codes — a one-dimensional version of the compressible code, 
CASTRO (JAlmgren et al.ll2010l ). in spherical coordinates; a one-dimensional version of MAESTRO 



-19- 



in spherical coordinates, and a full three-dimensional spherical star in MAESTRO. 

Our initial model is generated by specifying a core density (2.6 x 10^ g cm~^), temper- 
ature (6 X 10*^ K), and a uniform composition (X(^^C) = 0.3, X{^^0) = 0.7) and integrating 
the equation of hydrostatic equilibrium outward while constraining the specific entropy, s, 
to be constant. In discrete form, we solve: 

Po,j+i = Po,j + -Ar{poj + poj+i)g^+y^. (22) 

so,j+i = soj (23) 

We begin with a guess of Po j+i and Toj+i and use the equation of state and Newton- Raphson 
iterations to find the values that satisfy our system. Since this is a spherical, self-gravitating 
star, the gravitation acceleration, Qj+y^, is updated each iteration based on the current value 
of the density. Once the temperature falls below 10^ K, we keep the temperature constant, 
and continue determining the density via hydrostatic equilibrium. This uniquely determines 
the initial model. 

For the one-dimensional simulations, we map the inner 5 x 10® cm of the model onto 
a one-dimensional array with 1280 elements with Ar = 3.90625 x 10^ cm. For the full- 
star three-dimensional simulation, we map the model onto a 5 x 10® cm^ domain with 256^ 
Cartesian grid cells with Ax = 5Ar = 19.53125 x 10^ cm. For the one and three-dimensional 
MAESTRO calculations, we use cutoff densities (see Appendix I A. 5 1) of Pcutoff = 10^ g cm~^ and 
Paneiastic = 10^ g cm~'^. Corresponding to radii of approximately 1.8 x 10® cm and 1.9 x 10® cm, 
so the star easily fits within the computational domain for each problem. For the full-star 
three-dimensional simulation, we use an inner sponge with Vgp equal to the radius where 
Po = 10^ g cm"'^, Tind equal to the radius where po = 3 x 10^ g cm~^, and k = 10 s~^. We 
use the same outer sponge as in Paper IV. All boundary conditions are outflow, except for 
the center of the one-dimensional simulations, which uses a symmetry boundary condition. 
We run each simulation using a CFL number of 0.5. 

We heat the center of the star for 0.5 seconds and look at the solution at 2.0 seconds (after 
the compressible solution has had time to re-equilibrate). We use ifext = Hoe~^^^^^ '^^^ , 
with Hq = 10^^ erg g^^ s"^ chosen to be much larger than the nuclear energy generation 
rate during the convective phase of SN la, in order to see a measurable effect over a few 
seconds. A comparison of po,Po, and T at t = and t = 2 s for each code is shown in Figure 
[71 There is excellent agreement between each of the simulations. 
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6.3. Convergence Test 

In Paper III, we demonstrated that our single-level algorithm is second-order in space 
and time by tracking a hot bubble rising in a white-dwarf environment using two-dimensional 
planar geometry. Here, we perform the same test to show that our algorithm, with a level 
of refinement, is second-order in space and time. 

We choose a domain size of 7.2 x 10^ cm by 2.88 x 10^ cm and generate a high resolution 
initial model with Ar = 7.03125 x 10^ cm (which is equal to Ax for our highest resolution, 
"exact" solution) using the method described in Appendix O with rbase = 0. For each of the 
remaining, lower resolution simulations for this test, we generate an initial model using Ar 
equal to the effective Ax of each simulation by linearly interpolating values from the high 
resolution model. Next, we add a temperature perturbation of the form: 



Tij — Tq j + 0.3 



1 + tanh ( 2 - ^ 
a 



(24) 



where a = 2.5 x 10^ cm and dij is the physical distance between the cell center corresponding 
to cell {i,j) and the location (3.6 x 10^ cm, 3.2 x 10^ cm). Then, we call the equation of 
state to compute a consistent p,h = p, h{po, T, X^) everywhere. We use the reaction network 
described in § 4.2 in Paper III. We use cutoff densities of pcutoff = Panciastic = 3 x 10^ g cm~'^, 
and a sponge with Vgp equal to the radius where po = lOpcutoS, ''"md equal to the radius where 
Po = Pcutoff, and K = 10 s~^. We specify periodic boundary conditions on the side walls, 
outflow at the top, and a solid wall at the bottom of the domain. 

Since we do not have an exact analytical solution, we consider a single-level simulation 
run with 1024 x 4096 cells and At = 3.125 x 10"'^ s to be the exact solution. We perform 
three single-level simulations using resolutions of 64 x 256, 128 x 512, and 256 x 1024 grid cells 
using fixed time steps of At = 0.05 s, 0.025 s, and 0.0125 s, respectively. We also perform 
two simulations with a single level of refinement with effective resolutions of 128 x 512 
and 256 x 1024 grid cells with fixed time steps of At = 0.025 s and 0.0125 s, respectively. 
These fixed time steps correspond to a CFL of 0.9. We refine all cells in the range r G 
[1.8 X 10^, 5.4 X 10^] cm, so effectively we have refined 1/8"^ of the domain, making sure the 
hot spot is contained within the refined region. We run each simulation to t = 1 s. For this 
test, whenever we call Enforce HSE to compute po from po, we use r = 1.8 x 10^ cm as the 
starting point for integration rather than the location of the cutoff density to ensure that 
numerical errors due to integrating the equation of hydrostatic equilibrium across simulations 
with different resolutions are minimized. 

In order to compute the Li error norm for each simulation, we average the data from 
the exact solution onto a grid with corresponding resolution. We measure the Li error norm 
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in the physical space corresponding to the refined region using 

Ll = > \(f)ij — 0i7,cxact ! (25) 

where rZceii is the number of cells we sum over. This form of the Li error norm gives us 
the average error per cell. We compute the convergence rate, p, between a coarser and finer 
simulation using 

bl -^1, coarser \ /f->/^\ 

g2 -J • (26) 

\ -t^l,finer / 

Tables [1] and |2] show the Li error norms and convergence rates for the single- level and multi- 
level solutions, respectively. The convergence rates correspond to the two columns on either 
side of the reported value. We note second order convergence in each variable. Additionally, 
the magnitude of the Li error norms for the multi-level simulations is comparable to the 
corresponding resolution error norms for the single-level simulations. This means that the 
multi-level simulations are accurately capturing the finer-scale features, as compared to the 
single- level simulations, i.e., the presence of coarse grid data and/or coarse-fine interfaces is 
not harming the solution in the refined region. 



6.4. Adaptive Bubble Rise 

To test the ability for an adaptive, three-dimensional planar simulation to track a lo- 
calized feature, we examine a hot bubble rising in a white-dwarf environment. The problem 
setup is exactly the same as in §6.3[ except that we now compute in three dimensions and 
allow the grid structure to change with time. We choose a domain size of 7.2 x 10^ cm 
by 7.2 X 10^ cm 2.88 x 10^ cm and for each simulation, we generate an initial model with 



Table 1: Li error norms and convergence rates for the single- level simulations. 
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Ar = 5.625 x 10^ cm (which is equal to the effective Aa; for both of the simulations in 
this test) using the method described in Appendix O with rbase = 0. We add a tempera- 
ture perturbation of the form given in equation fl2^ . but in three dimensions with the hot 
spot centered at location (3.6 x 10^ cm, 3.6 x 10^ cm, 3.2 x 10^ cm). We will show that 
the adaptive simulation captures the same dynamics as the single-level simulation in a more 
computationally efficient manner. 

We compute a single-level simulation with 128 x 128 x 512 grid cells, and an adaptive 
simulation with two levels of refinement at the same effective resolution. For each cell, if the 
T — T > 3 X 10^ K, we tag all cells at that height to ensure they are computed at the finest 
refinement level. We run each simulation to t = 2.5 s using a CFL number of 0.9. 

Figure [H] shows the initial profile oiT — T and the initial grid structure of the multi- level 
run. The single-level simulation has 8,388,608 grid cells and takes approximately 32 seconds 
per time step. The adaptive simulation initially has 131,072 grid cells at the coarsest level, 
114,688 cells at the first level of refinement, and 688,128 cells at the finest level of refinement 
(the number of grid cells at the finer levels changes slightly with time as the grid structure 
changes to track the bubble) and takes approximately 12 seconds per time step, for a factor 
of 2.7 speedup. Both computations were performed using 32 processors on the Franklin XT4 
machine at NERSC. Figure [9] shows a series of planar slices of the simulations at 0.5 second 
intervals in order to show that the adaptive simulation captures the same dynamics as the 
single-level simulation. The vertical distance shown is from z = to 9.2 x 10^ cm. 

Table 2: Li error norms and convergence rates for the multi-level simulations. 
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6.5. Forced Convection 

To test the expansion of the base state in a multi-level, two-dimensional planar sim- 
ulation, we simulate a white-dwarf environment with an external heating layer. We also 
show that the multi-level simulation captures the same fine-scale structure as a single-level 
simulation at the same effective resolution for a short time. We performed a similar test in 
Paper III, but without refinement. 

We choose a domain size of 2.5 x 10^ cm by 4 x 10^ cm and for each simulation, we 
generate an initial model with Ar equal to the effective Ax for that simulation using the 
method described in Appendix O with rtase = 5 x 10^ cm. The low entropy layer beneath 
our model atmosphere prevents the convective flow from reaching our lower boundary. 



We apply an external heating layer of the form 



i7ext,«- = i/oe^'^^-'^'--)' 



m=l \ ^x / 



(27) 



with riaycr = 1-25 x 10^ cm, Hq = 2.5 x 10^^ erg g ^ s ^, L^ = 2.5 x 10® cm, and Vj and Xi being 
the radial and horizontal physical coordinates of cell {i,j)- The perturbation parameters are 
b = (0.00625,0.01875,0.0125), A; = (2,6,8), and ^ = (0, 7r/3, 7r/5). We disable reactions for 
this test, since the heating layer was chosen to expand the base state over a very short time 
period, rather than accurately model reactions. 

We use cutoff densities of Pcutos = Panciastic = 3 x 10^ g cm~^ and a sponge with rgp 
equal to the radius where po = 10® g cm~^, r^d equal to the radius where po = Pcutofr, and 
K = 100 s~^. We specify periodic boundary conditions on the side walls, outflow at the top, 
and a solid wall at the bottom of the domain. 

In this test, we use a CFL number of 0.9. We perform two single-level simulations using 
80 X 128 and 320 x 512 cells, and a simulation with two levels of refinement and an effective 
resolution of 320 x 512 cells. For the multi-level simulation, we fix the refined grids, ensuring 
that r G [9.375 x 10^, 1.5626 x 10®] cm (which contains the external heating layer) is at the 
finest level of refinement. We run each simulation to t = 30 s. 

Figure [TO] shows po,Po, and T after 30 seconds of convection for the single-level fine grid 
simulation and the multi-level simulation. There is an excellent agreement between these 
two simulations, except in the temperature profiles at the top of the domain. However, 
this corresponds to a region with density below the cutoff densities where the temperature 
is extremely sensitive to small density perturbations, and furthermore, is not fully refined 
in this test, so this is an acceptable difference. Both simulations were performed using 4 
processors; the single-level fine grid run required approximately 1.9 seconds per time step 
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and the multi-level run required approximately 0.6 seconds per time step, for a factor of 3 
speedup. 

Figure [TT] shows the temperature profile after 3 and 4 seconds for each of the three 
simulations. The vertical distance shown is from z = 5 x 10^ cm to 2.2 x 10^ cm. At t = 3 s, 
the multi-level simulation is able to capture the finer-scale structure visible in the single- 
level fine grid simulation, which is not captured in the single-level coarse grid simulation. At 
t = 4 s, a finer-scale structure is still visible in the multi-level simulation, but the solution 
begins to drift from the single-level fine grid simulation, which is expected since we are 
deliberately refining only a part of the convective region. 



6.6. Full-Star AMR 

We now compute three-dimensional, full-star calculations of convection in a white dwarf. 
This problem models the convection and energetics of a white dwarf that is a few hours from 
reaching ignition. We performed similar simulations using an earlier version of the algorithm 
in Paper IV. 

We begin with the one- dimensional white dwarf model described in §2.4 in Paper IV. We 
map this one-dimensional model into the center of a computational domain of 5 x 10*^ cm on 
a side. The first simulation is a single-level simulation with 384'^ cells. The second simulation 
is adaptive with two levels of refin ement and an eff e ctive resolution of 384^ cells. We use 



the reaction network strategy from IChamulak et al.l (120081 ) to compute the energetics from 



the carbon burning. This modified network differs from the one used in paper IV in the ash 
composition (we now burn to an ash with A = 18 and Z = 8.8) and the energy re lease (we 



use a quadratic fit to the g- values tabulated on page 164 of IChamulak et al.ll2008l ). Finally 
instead of destroying two carbon nuclei for each reaction, we use the Mu value of 2.93 
describe in that paper. We initialize the simulation with a velocity perturbation described 
exactly as in §2.4 in Paper IV. We use the same cutoff densities, sponge parameters, and 
boundary conditions as in §6.2[ 

Figure [12] shows the initial grid structure of the adaptive simulation. Based on our work 
in Paper IV, we choose to fully refine all cells where p > 10® g cm~'^, since at early times, 
the dynamics of the star are driven by the reactions and convection in this inner region. We 
wish to examine whether the adaptive simulation can give the same result as the single-level 
simulation, and the computational efficiency of each simulation. 

We use a CFL number of 0.7 and compute to t = 900 s. We choose two diagnostics used 
in Paper IV to compare the simulations. Peak temperature is a useful diagnostic since the 
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reaction rates are extremely sensitive to temperature, and thus peak temperature serves as 
a good guide for observing the progression toward ignition. Peak radial velocity is another 
useful diagnostic as it is a simple measure of the strength of the convection within the star. 
Since the solution of ou r equation is highly non-linear, and the reaction rates scale with 



temperature as ~ T^^ (jWoosley et al.ll2004l ). we expect that errors from the coarse grid 



will perturb the solution at the finest level, eventually causing significant differences in the 
exact flow field. However, as shown in Paper IV, when we run our simulation to the point 
of ignition, we require upwards of hundreds of thousands of time steps. Therefore, in the 
comparison diagnostics in this test, it is sufficient to compare the overall qualitative solution. 
An exact quantitative comparison is impossible over long times. 

Figure US] shows the evolution of the peak temperature and peak radial velocity over 
the first 900 s for both simulations. The adaptive simulation gives the same qualitative 
result as the single-level simulation. As mentioned before, the curves do not match up more 
closely because the equations are highly non-linear, and slight differences in the solution 
caused by solver tolerance and discretization error change details of the results, but not the 
overall qualitative results. The single-level simulation has 56,623,760 grid cells and takes 
approximately 36 seconds per time step. The adaptive simulation initially has 884,736 grid 
cells at the coarsest level, 3,511,808 cells at the first level of refinement, and 4,282,048 cells 
at the finest level of refinement (the number of grid cells at the finer levels changes slightly 
with time) and takes approximately 18 seconds per time step, for a factor of 2 speedup. Also 
note that the overall memory requirements are significantly less for the adaptive simulation, 
as can be seen by the reduced total cell count. Each simulation was run using the Franklin 
XT4 machine at NERSC with 216 processors. 

We note that in the future, when we perform longer calculations up to the point of 
ignition, we may have to refine a greater portion of the star in order to properly capture 
the overall dynamics as the convective region expands. In a forthcoming paper, we plan on 
performing higher resolution studies, where AMR will save us both time and computational 
resources. 



7. Conclusions and Future Work 

We have developed a low Mach number hydrodynamics algorithm suitable for full star 
fiows and local atmospheric regions with a time-evolving base state within an AMR frame- 
work. In forthcoming papers, we will use MAESTRO to further our scientific investigation 
of the convective phase of SNe la. Our previous simulations in Paper IV were at modest 
resolution and assumed a constant base state. We are now performing simulations at higher 
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effective resolutions with the use of AMR along with a time-varying base state. As part of 
this study, we are examining the tagging conditions necessary to model a full-star up to the 



point of ignitio n. We are also studying Type I X-ray bursts fIStrohmayer fc Bildstenl 12006 



Lin et al.ll2006l ). which are believed to be caused by the thermonuclear explosive burning of 
hydrogen and/or helium gas accreted into a thin shell on the surface of neutron stars. We 
pose this problem in planar geometry, model a small patch of the star, and refine near the 
base of the accreted layer. 
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A. Time Advancement Algorithm 

A.l. Single-Level Algorithm Changes From Papers III and IV 

The single-level algorithm has gone through numerous changes since Papers III and IV. 
The current single-level algorithm is presented in §A.2f[X75| we summarize here the changes 
since Papers III and IV: 



We have extended the base state evolution to spherical problems by defining Advect 
Base Density Spherical ( §A.3.2I) . Advect Base Enthalpy Spherical ( §A.3.4p . 
and Compute wq Spherical ( §A.3.5|) . We have also defined spherical discretizations 
for ^ and r]p (Steps 4C, 4F, 8C and 8F in ^ OOj) . 

For spherical problems, we have improved the accuracy of the mapping of data between 
the one-dimensional radial array and the three-dimensional Cartesian grid (SD. 

We have upgraded our unsplit piecewise-linear Godunov method for time-centered edge 
state prediction for the ba se state and Cartesian grid data to use the unsplit piecewise- 
parabolic m ethod fPFM) fjCq 



dimensions (JMiller fc Colella 



ella fc Woodwardlll984f) with full corner coupling in three 



2OO2I : ISaltzmanlll994l ). We shall refer to this procedure 
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as "computing the time-centered edge states" ( §A.3.2t §A.3.4t and Steps 3, 4, 7, 8, 
and 11 in ^KM. 



As introduced in Paper IV, we update T after the call to React State ( §A.3.ip . 



Rather than evolve po using an evolution equation, we simply update po using the 



hydrostatic equilibrium equation ( §A.3.3p . These two methods are analytically equiv- 
alent, but in our experience, the numerical drift associated with evolving po using an 
evolution equation causes the entire solution to drift from thermodynamic equilibrium 
over time. 

As explained in ^ in the advection step, we predict (ph)' to edges instead of T 
(Steps 4Hi and 8Hi in §A.4p . Thus, a base state enthalpy, {ph)o is required in order 
to construct the enthalpy fluxes for the conservative update. Unlike po, we do not 
need to carry the complete evolution of (p/i)o- In practice, we set {ph)o equal to 
the lateral average of the full state enthalpy after the first call to React State, i.e., 
(p/i)o = (p/i)(^). We then advect {ph)^ without reactions or heating to mirror the 
treatment of (ph) in the advection step (Steps 4G and 8G in §A.4I) . 

In the advection step, rather than simultaneously updating each of the base state 
quantities, followed by an update of all the full state quantities, we have interwoven 
these updates in order to obtain better accuracy. In order: we advect po, advect p, 
correct pq, advance po, advect {ph)o, and advect (ph) (Steps 4 and 8 in §A.4p . This 
enables us to use a time-centered ip rather than a time-lagged tp for the enthalpy 
update. 

Rather than compute V-?7p and use it to adjust po to account for convecting overturning, 
as was done in Correct Base in Paper III, we simply use the lateral average operator 
to enforce po = P, which is simpler and analytically equivalent (Steps 4D and 8D 
in §A.4[) . We use the improved averaging procedure described in §4.11 which greatly 
improves the accuracy of this mapping for spherical problems. 

We have moved the first reaction step (formerly Step 3 in §A.4I) to occur before Steps 
1 and 2 in §A.4[ This is a non-functional change; it was only made so that Steps 2-5 
mirror Steps 6-9 in the overall predictor-corrector scheme. 

For planar problems, we evolve the enthalpy for the sole purpose of updating the tem- 
perature, which is subsequently fed into the reaction network and used in computing 
thermodynamic derivatives. For spherical problems, as described in Paper IV, we in- 
stead define T from p,po, and Xk. This completely decouples the enthalpy from the 
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rest of the solution. Our experience has shown that, with spherical geometry, the dis- 
cretization errors are minimized by using the hydrostatic, radial base state pressure 
to define temperature. We still evolve the enthalpy, but we use it as a diagnostic to 
examine to what extent our numerical method keeps the solution in thermodynamic 
equilibrium. 



A. 2. Gravity 

For planar problems, we treat gravity as constant in space and time. For spherical 
problems, gravity is computed directly from pQ in the following manner. First, we define the 
enclosed mass at radial edges and cell centers: 

fc=l 

In practice, we compute (r^_i, —rl_3, ) as Ar(r^_i, +rf.__y^rf._3/^ + r^_a, ), in order to minimize 
roundoff error at large radii and use a similar formula for the radial difference term in the 
equation for rricncij- Next, we define the gravity at both radial edges and cell centers: 

93-% = ——2 ; 9 J = ■ (A2) 

We indicate which base state density is used to compute gravity by using a shorthand 
notation with superscripts on g, e.g., g"" = gipo)- 



A. 3. Main Algorithm Notation 

We make use of the following shorthand notations in outlining the algorithm: 

A.3.1. Reactions 



React State[p'°, {ph)"^ , X]^ ,T'^ , (pffext)'",p|,1 ^ [p°"S (p/i)°"*,Xr*,T°"\ (pcO^)""*, {pH, 
evolves the species and enthalpy due to reactions through At/2 according to: 



\outl 



nuc/ 



dXk . , . dT 

——=Uk{p,Xk,Ty, -— 

at at c, 



-{-Y^^kCok + H^ncY (A3) 
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Full details of the solution procedure can be found in Paper III. We then define: 



{pCokT' 



{phf 



At/2 



°- + ^(pi7ext) 



(A4) 
(A5) 



where the enthalpy update includes external heat sources (pi/cxt)™- As introduced in Paper 
IV, we update the temperature using T°"* = T{p°'^^, h°'^^, X^^^) for planar geometry or T°"* = 
T(p°"*,Pq'^,X^"*) for spherical geometry. Note that the density remains unchanged within 
React State, i.e., p°"* = p'°. 



A. 3. 2. po Advection 



Advect Base Density[p|f , wj,^] -^ [p§"*,Po 



,out out,n+y2-| . 



is the process by which we update 



the base state density through At in time. We keep the time-centered edge states, Pq^ '" , 
since they are used later in discretization of r]p for planar problems. 



planar: We discretize equation fll4p to compute the new base state density, 



out _ in 
P0,j — HO,j 



At 
Ar 



out,n.+y2 in 



Po 



w, 



j+'h 



out,?i+y2 



Po 



Wn 



i-% 



(A6) 



We compute the time-centered edge states, Pq" '" , by discretizing an expanded form 
of equation ([1 



'9po 9po 

dt dr 



-po- 



dr 



(A7) 



where the right hand side is used as the force term. 
spherical: The base state density update now includes the area factors in the divergences: 



out _ in _ ]_^ 
Po,j - Po,j ^2 ^^ 



^2^out,n4-y2^in 



i+y2 



^2^out,n+y2^in 



J-% 



{Al 



In order to compute the time-centered edge states, an additional geometric term is 
added to the forcing, due to the spherical discretization of ( 1T^ : 



dpo dpo 



dt 



Wo- 



dr 



-po- 



dvjQ _ 2poWo 
dr r 



(A9) 
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A.3.3. po Update 



Enforce HSE[p™,p™] — )■ [po"*] has replaced Advect Base Pressure from Paper III 
as the process by which we update the base state pressure. Rather than discretizing the 
evolution equation for pq, we enforce hydrostatic equilibrium directly, which is numerically 
simpler and analytically equivalent. We first set Po"/=o ~ Po'i=o ^^^ then update pg"* using: 



Po">i = pZ + ^^^.+% 



fc+i + p^:,) 



(AlO) 



where g = g{pQ^). As soon as p^j < PcutoS, we set Poj'+i — Po^j ^^^ ^^ remaining values of j. 
Then we compare Pq^}^ with p'i^j^^^^^ and offset every element in pg"* so that Po"imax ~ Po'jmax' 
We are effectively using the location where the p™ drops below Pcutoff as the starting point 
for integration. 



A. 3. 4- iph)o Advection 

Advect Base Enthalpy[(p/i)™, Wg^-?/;™] — > [(p/i)o"*] is the process by which we update 
the base state enthalpy through At in time. 



planar: We discretize equation (fT5|) . neglecting reaction source terms, to compute the new 
base state enthalpy. 

At 



(Phr- = (Ph) 



0,3 



Ar 



{ph)TM 



3 + 'h 



{phT.^M 



J i-y2 



Atiijf. (Aii: 



We compute the time-centered edge states, (p/i)o , by discretizing an expanded form 
of equation (fTSl) : 



d{ph)o d{ph)o 



-iph)o^r- + 'ip- 



(A12) 



dt dr dr 

spherical: The base state enthalpy update now includes the area factors in the divergences: 



{phX- = iph)Z 
1 At 



r? Ar 



'r\phC'/^w-^^ 



Jj+y2 



r\ph)r'< 



J j~Y2 



At^Jj 



in,n+y2 



(A13) 



In order to compute the time-centered edge states, an additional geometric term is 
added to the forcing, due to the spherical discretization of ( IT5l) : 



(AM) 



d{ph)o , d{ph)o dwo 2{ph)oWo 



dt 



dr 



dr 
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A. 3. 5. Computing Wo 

Here we describe the process by which we compute wq- The arguments are different for 
planar and spherical geometries. 

Compute Wo Planar [5™, r\''',p|jW'1 ^ K^*]: 

In Paper III, we showed that ip = rjpQ for planar geometries, and derived an alternate 
expression for equation (fTTI) . We compute Wq using equation (35) in Paper III using 
the following discretization: 

with WQ_y = 0. 



Compute Wo Spherical [S , Fi ,p^,p^,r]'^] -^ K"*]: 



We begin with equation ( uTj) written in spherical coordinates: 



(.^,„.,„).,„(5__i_|^). (A16) 



We expand the spatial derivative and recall from Paper I that 

1 dpo _ 1 d^Q 
I\po dr /3o dr ' 

giving: 

1 «9 / 2 \ -n 1 f9po . dpo 



(A17) 



r^dr^^'-^^-^~Wo[^^^'^^- ^^^'^ 



We solve this equation for wq as described in Appendix | 



A. 4. Main Algorithm Description 



We now describe the main algorithm, making frequent use of the shorthand developed 
above. In summary, in the predictor step (Steps 2-5) we use an estimate of the expansion 
term, S, to compute a preliminary solution at the new time level, denoted with an "?7, + 1,*" 
superscript. In the corrector step (Steps 6-9), we use the results from the predictor step to 
compute a more accurate expansion term, and compute the final solution at the new time 
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level, denoted with an "n + 1" superscript. We use Strang-splitting to achieve second-order 
accuracy in time. See Figure [Hlfor a flow chart of the algorithm, including the notation used 
as we advance the solution by At. Figure [T^ is a flow chart of the advection steps (Steps 4 
and 8), which includes the notation we use as we advect the solution through a time interval 
of At. 

The discussion that follows mirrors closely that in Paper III, but has been updated to 
reflect all the changes throughout the algorithm. The advance of the state through a single 
time step appears as: 

Step 0. Initialization 

This step remains unchanged from Paper III. The initialization step only occurs at the 
beginning of the simulation. The initial values for \J^,p^, {ph)^,X^,T^, Pq,Pq, and P'j' 
are specified from the problem-dependent initial conditions. The initial time step, At°, 
is computed as in Paper III. Finally, initial values for w^^ , rjp , 4>~^'^, 7r~/^ 5'°, and S^ 
come from a preliminary pass through the algorithm. 

Step 1. React the full state through the first time interval of At/2. 

Call React State[p", (p/i)",X^,T", (p/7ext)",Po] 

-> [p«, (p/.)«,xW,TW, (pa;,)W, ipH^a% 



Step 2. Compute the provisional time-centered expansion, 5'""^/^'**, provisional base state 
velocity, Wq , and provisional base state velocity forcing. 

A. Compute S'^'^^'^'**. We compute an estimate for the time-centered expansion term 
in the velocity divergence constraint (eq. [I2]). For the first time step {n = 0), 
we set 

where S^ i s foun d during initialization. For other time steps (n 7^ 0), following 



Bell et al.l ( 120041 ). we extrapolate to the half-time using S at the previous and 

current time levels 

At"" ^" - ^"-1 

At^ 



^n+%.. ^gn^ f^:l__i_. (A20) 



Next, compute 



B. Compute Wq . 



Sn+%^^ = Avg (s''+'b'**\ . (A21) 



For planar geometry, call 

Compute Wo PlanarfS^+y^.*'^, If, p[}, ^"-V^] -^ K'^'''''^] 



-33- 



For spherical geometry, call 



Compute Wo Spherical[5"+y2'**, r^,p^,p^, r/p ''] -^ [w, 



"■-y^] ^ r^,,"+y2,*i 







C. Compute the provisional base state velocity forcing. Rearrange equation ([9]), 

1 dTTo dwo dwo , . 

+ «^o^r-, (A22) 



Po dr dt dr 

with the following discretization: 

^^^°^ ^0 -^0 <-f^) , (A23) 



^po dr J (At" + At"-i)/2 " \ dr 

where Wq'* and (dwo/dr)^'* are defined as 

<•* = ^^^° +^/ f° , (A24) 

At" + At"-i ^ ^ 






A,.|^V-\a,.-.P»o^"^'"* 



9r y At" + At"-i " \ dr J ' " \ dr 

(A25) 

If n = 0, we use At"^ = At°. 

Step 3. Construct the provisional time- centered advective velocity on edges, XJ'^^^'*. 

Using equation (fTOl) . we compute time-centered edge velocities, IJ'^^^'"'''*, using U = 
U" + Wq . The t superscript refers to the fact that the predicted velocity field does 
not satisfy the divergence constraint. We then construct XJ^°^'* from XJ^^^'^'* using 
a MAC projection, as described in detail in Appendix B of Paper III. We note that 
-y-ADV,* ga^igggg w^Q discrete version of (U"^°^'* ■ e^) = as well as 

V ■ U^V^^^^A = /3o" fs''+'^^'** - S^'+'k'**) , (A26) 

/3o = /3o(Po,Po,rT), (A27) 

where /3o is computed as described in Appendix C of Paper III. 
Step 4. Advect the base state and full state through a time interval of At. 

A. Update po, saving the time-centered density at radial edges by calling 
Advect Base Density [p^, wio"""^"'] ^ [p?"^'^P^'''"*•"''1. 

B. Update (pX^) using a discretized version of equation ([1]) omitting the reaction 
terms, which were already accounted for in React State. The update consists of 
two steps: 
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i. Compute the time-centered species edge states, (pXfc)"'+/2:*.pred^ ^^^ ^^^q con- 
servative update of {pXi^Y^\ We use equations flT7|) and flT3|) to predict 



p\-/ = p{i) _ p^ g^T^fi x^' = (pXf^Y^yp^^^ to time-centered edges using 



V(i) 



w, 



n+%,* 



e„ yielding p'«+y2:*>pred ^nd x^+'/^'*'P'''''^. We convert 



the perturbational density full state density using 



n+y2,*,prcd _ 'n+y2,*,prcd , Po + Po 



(2a), 



(A28) 



where the base state density terms are mapped to Cartesian edges. Then, 

ioX ')"+y2.*,pred _ / n+y2,*,prcd-v-"+y2i*.PrcdN 



ii. Evolve (pX/, 



^(1) 



(P^A 



U2),* 



usmg 



{pX.t'^^^ 






U^°^'* + t.o"^^/-V)(pX,)' 



\n+y2,*,prcd 



P 



(2),^ 



5^(px,)(2)'^ xf '^ = (px,)(^)'7p(^)'^ 



(A29) 
(A30) 



C. Define a radial edge-centered rfl, '^'*. 



For planar geometry, since rfp = p'(U ■ e^) = p(U ■ e,.) — po(U ■ e^.) = p(U ■ e^/ 
PoM^o, 



^n.+y2,* 



Avg$:f(u 



tADV.* 



e^ + w, 



n+y2,*\ /'^ V. ^"•+y2,*,prcd 







(pXfc 



„, "+y2,* ^i+y2,*,prGd 

"^0 Po ' 



(A31) 



For spherical geometry, first construct rfi^^ '" = [p'(U-er)]"'^/2'* on Cartesian 
cell centers using: 



cart,n+y2,* 
Ip 



pW + p(2),^ 



2 

TjADV,Vr 



Po+Po 






Then, 



^n+y2,* _ ^^g j'^cart,n+y2,*j _ 

ntered ?7p . To get 
the two neighboring radial cell-centered values 



(A32) 
(A33) 



This gives a radial cell-centered ?7p . To get ?7p at radial edges, average 
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n+l,*^ 



D. Correct po by setting Pq"*" '* = Avg(p'^^)'*). 

E. Update po by calling Enforce HSE[p[J,po+^'''] -^ [p^ j. 

F. Compute ^"+'/2'*. 

For planar geometry, 

For spherical geometry, first compute: 



ri 



(1) 



r\ 



(2),^ 



Avg 
Avg 



J- 1 \P SPO'^fc 



r, p(^)-,pr^^4^ 



(2),* 



(A34) 

(A35) 
(A36) 



Then, define ip'^+z^'* using equation (lAlSp 



^ 



n+y2,* 



rji) + r(2). 



pS + Po 



5: 



n+y2,* 



2^,,r!-+y2,* 



r W, 







i+y2 



2 "+y2,- 



j-'h 



(A37) 



G. Update {ph)o. First, compute {ph)^ = Avg[(p/i)(^^]. Then, call 
Advect Base EnthaIpy[(p/i)^,M;o+'/''*,7/^"+'/2,^] ^ [(p/i);;+i'*]. 

H. Update the enthalpy using a discretized version of equation (J3]), again omitting 
the reaction and heating terms since we already accounted for them in React 
State. This equation takes the form: 



djph) 
dt 



-V-(Up/i)+^ + (U-e,)^. 

or 



For spherical geometry, we solve the analytically equivalent form, 
d{ph) 



dt 



-V ■ (Vph) +^ + V- (Vpo) - PoV ■ U, 



(A38) 



(A39) 



which experience has shown to minimize the drift from thermodynamic equilib- 
rium. The update consists of two steps: 

i. Compute the time-centered enthalpy edge state, (p/i)"-+/2'*'Pi''=d^ for the con- 
servative update of {phY^\ We use equation (fT8|) to predict (ph)' = (ph)^^^ — 
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(p/z)o to time-centered edges, using U = \J^^^'*+WQ~^'^'*er, yielding (p/i)'"+y2>*,prcd_ 
We convert the perturbational enthalpy to a full state enthalpy using 

For planar geometry, we map (p/i)o directly to Cartesian edges. In spherical 
geometry, our experience has shown that a slightly different approach leads to 
reduced discretization errors. We first map ho = (p/i)o/po and po to Cartesian 
edges separately, and then multiply these terms to get (p/i)o- 
ii. Evolve (p/i)(i) -^ (p/^)(2).^ 
For planar geometry, 

-At I V ■ [(U^°^'* + u;o"^'/"*e,) (p/i)"+y^'*'P^^'i] } 
+At (U^°^'* ■ e,) (^^ " + At^"+y-^ (A41) 

For spherical geometry, 

-At {v • [(u^°V'^ + u'o+'/^'^e,) (p/i)"+y^'*'P^'='^] } 

+At { V ■ (U^°^'X) - Po V ■ U^°V'*} 

+AtV'"+'/^•^ (A42) 

Then, for each Cartesian cell where p^^)-* < pcutoff, we recompute enthalpy using 

(p/.)(^)'* = p(^)'^/.(p(^)^pr^^xf'^). (A43) 

I. Update the temperature using the equation of state: T'^^)'* = T{p^'^^'*, /i'^^^'*, X^ ) 
(planar geometry) or T'^^^'* = T(p'^2)'*,pQ '*,X]. *) (spherical geometry). 

Step 5. React the full state through a second time interval of At/2. 
Call React State[p(2).^ (p/i)(2).*,xf •*,T(2).*, {pH,^tY^^'\p^^''*] 

^ [p"+l•^ (p/^)"+l.^ x,"+l'^ ^"+l.^ (pcI;,)(')'^ (pi/nuc)(')'i. 

Step 6. Compute the time-centered expansion, 5*"+/^'*, &ase state velocity, Wq , and base 
state velocity forcing. 
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A. Compute 5*"+/^'*. First, compute S"'~^^'* with 

^n+i. = -aJ2U^>^r^' + -^ Y.P^.M^"^''' + ^H^ + -^S^ (A44) 
k P ' Pp k 

where (cufc)*^^^'* = (pcok)^'^^'* /p^'^^'* and the thermodynamic quantities are defined 
using p"+^'*, X^~^ '*, and T"+^'* as inputs to the equation of state. Then, define 

en I cn+l,* 

^n+%,. = Avg(5"+'/2'*), s'^+'Z'-'' = ^—^ , (A45) 

B. Compute Wq . First, define 



r"+/2.* _ -*- 1 + ^ 1 ^n+Ya,^ _ Po + Po ^rt+Va,* _ Pq + Pq I \ ar\ 

^ 1 — 7, ' Po — 7, ' Po — 7 ' l,A40j 



with 



rr+i'* = Avg [Fi (p"+^'^p^''^x,"+l'*)] . (A47) 

For planar geometry, call 

Compute Wo V\axvBT\S^+^Xi^''^'\pV'''''\V^''''*] ^ V%^H 
For spherical geometry, call 

Compute Wo Spherical [:S^^^%^, ^^'/^'^ p^'/^'^p^'/^'^ r/r'^"*] 

r n+y2T 

^K ]• 

C. Compute the base state velocity forcing. Rearrange equation (]A22p . 



ldn,Y_ wT'^^-wr''^ w-J^-^] , (A48) 



^po dr J i/2(At" + At"-i) ° V dr 

where Wq and (dwo/dr)"' are defined as 






^ ^^^0 +At i/^o (A49) 

" At" + At"-i ^ ^ 



9r / At" + At'^-i 



Ar(^V"' + Ar-^^^°^'^^' 



9r / \ dr 



(A50) 
If ra = 0, we use At"^ = At°. 
Step 7. Construct the time- centered advective velocity on edges, U^^'^. 
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The procedure to construct U^^^'^ is identical to the procedure for computing U^°^'^'* 
in Step 3, but uses the updated values Wq and tTq rather than Wq and tTq'*. We 
note that IJ'^^^ satisfies the discrete version of (U'^^^ ■ e^) = as well as 

V ■ (/3o"+y^'*U^°v) = C^^"- [S^^'l^^^ - ;S^^^^) , (A51) 

on _|_ on+l,* , , 

Po 2 ' ^0 -PoyPo ,Po ,ii J- [AbZ) 

Step 8. Advect the base state and full state through a time interval of At. 

A. Update po, saving the time-centered density at radial edges by calling 
Advect Base Density [p-, 1/;^+^^] ^ [p^o"\ Po'''^"'^'"% 

B. Update (pX^). This step is identical to Step 4B except we use the updated values 
Wq , U^^^, and Pq rather than Wq , XJ^^^'*, and pg '^ . In particular: 

i. Compute the time-centered species edge states, (pXfc)"'"''/2:pred^ ^^^ ^j-^g ^q^_ 
servative update of {pXkY^\ We use equations flT7|) and flT3|) to predict 
p'(i) = pW — p^ and Xf. = (pXk)^^^ /p^^^ to time-centered edges with U = 
UADV ^^n+y2g^^ yielding p^+'h'P'^'^ and x^^'^^'"^""^ . We convert the perturba- 
tional density to a full state density using 

n _|_ (2a) 
n+i/a.prcd ^ 'n+Va.prcd _^ Po + Po (A53) 



Then, (pXfc)"+y2,prcd _ (^n+y2,prcdj^^"+y2,pred^_ 

ii. Evolve (pXfc)(i) -^ (pXfe)(2) using 

(pX,)(2) = ipX.p - At { V ■ [(U^°^ + ^o"^^''^e,) ipX.r^'/-^-'] } , (A54) 

P^^^ = E(P^^)^'^ ^^'^ = (p^^O^^Vp^^^- (A55) 

k 

C. Define a radial edge-centered r]^ . 
For planar geometry, 

<+y^ = Avg Y. [(U^°v ■ e,. + wT'l'^) {pX,r^'/-^-'' 

k 

-wI-^'^'pI-^"''^^''", (A56) 
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For spherical geometry, first construct r/p^"" '" = [p'(U ■ ef,)]"^/^ on Cartesian 



cell centers using: 



cart, 71+1/2 
'P 



p(l) + p(2) 



Po + Po 



(2a)' 



u 



ADV 



Then, 



7;+y^ = Avg ( 



cart,n+y2 
Ip 



e^ + wio'+'^n. (A57) 



(A58) 



This gives a radial cell-centered ?7p . To get r]^ at radial edges, average 
the two neighboring cell-centered values. 

D. Correct po by setting Pq^^ = Avg(p(^)). 

E. Update po by calling Enforce HSEK,p;j+i] -^ [p'^+'^]. 

F. Compute ^^"+^2. 

For planar geometry. 



1 ^^«+% 
For spherical geometry, first compute: 



r" 



2 l^p-i-y^ + \A) 9- 



Then, define ?/'""'" '^ using equation (1A18I) : 

-n+l/s _ 



(A59) 



(A60) 



^(1) 



^ 



pi^; + ri' 



(2)' 



r.n+y2 -■- 



Po+Po 



2„, "+y2 



n+l 



r w. 







2„, "+y2 



r U7, 







^3+y2 V /j-y2. 

G. Update (p/i)o by calling Advect Base Enthalpy [(p/i)[f, Wq^'^\ ip''+^ -^ [(p/i)o^^]. 

H. Update the enthalpy. This step is identical to Step 4H except we use the updated 
values Wo+'/^UA°v^^n+l^ (p/^)[5+^p[;+'/^ and ^P^'+'Z^ rather than 
^^+'k'\UADY,*^pn+i,.^ (p/i);;+i'*,p^, and tlj^+'/^'\ In particular: 

i. Compute the time-centered enthalpy edge state, (ph)^'^^'^'^^'^'^, for the conser- 
vative update of {ph)^^\ We use equation ( ITSl) to predict (ph)' = (ph)^^^ — 
(ph)^ to time-centered edges with U = IJ^^^+w'^^'^^Br, yielding {ph)'''+'/^'P''"^. 
We convert the perturbational enthalpy to a full state enthalpy using 



iph) 



n+y2,prcd 



iph) 



'n+y2,prcd , (P^)o + {ph)o 



n+1 



(A62) 
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ii. Evolve (p/i)(^) -^ {phY^\ 
For planar geometry, 



(p/i)(2) = (p/i)«-At{v-[(u^°^ + i/;o""'V)(p/^) 



n+V; 



'2 



tAdv .^\ I^\ , AfW,-+y2 



+At ( U^^^ ■ er ) ( ^ ) + At^'^+/^ (A63) 



For spherical geometry, 

{phY'^ = {phf^ - At {V ■ [(U^°^ + i/;o'"'^'e,) (p/i)"^^-^^^'^] } 



+At [V ■ (U^°VJ+/^ j - p^/^V ■ U^°^J + Atr^'/\ 

(A64) 

where Po^'^' is defined as p^^'^' = {p^ + p]^+^)/2. 
Then, for each Cartesian cell where p*^^-* < Pcutoff, we recompute enthalpy using 

(p/.)(^) = p(^)/.(p(^),pr\xf). (A65) 

I. Update the temperature using the equation of state: T*^^) = T{p'^'^\h'^'^\x\^ ) 
(planar geometry) or T*^^) = T {p^'^\ Pq'^^ , X\, ') (spherical geometry). 

Step 9. React the full state through a second time interval o/ At/2. 
Call React State[p(2), {phY^\ xf\T^^\ {pH,,tY^\p^+'] 



— > 



[p«+i, (p/^)«+i,Xr\T"+\ {pu,Y'\ ipH, 



nuci 



1(2) 



Step 10. Define the new time expansion, 5"""^^, and F'""''^. 
A. Define 

S-+' = -aY,U^kY'^ + ctH^^c + -^ Ep^^(^'^)^'^ + ^^e2, (A66) 
k P Pp k 

where {cOkY"^^ = (pi^fc)''^Vp*^^'* ^^^ ^^e thermodynamic quantities are defined using 
pn+i^ j^n+1^ ^^^ rpn+1 g^g inputs to the equation of state. Then, compute 



Sn+i = Avg(5"+^). (A67) 

B. Define 

F^i = Avg [Fi (p"+\pr\xr')] . (A68) 



-41- 



Step 11. Update the velocity. 

First, we compute the time-centered edge velocities, \j'^+ h^P'^^'^ ^ Then, we define 

P-^''^ = '-^. pT'^ = '-^^. (A69) 

We update the velocity field U" to U"^^'''" by discretizing equation ( ITO!) as 



U"+i.t = U"-At 



(^U^DV + w^+'ker'^ ■ VU^+'/^'P'^'^^' 



or 



dwr 



n+y, 



'2 



+At 



i.o.^-nfi^re,.- ^-"".:f") ..... 



pn+y2 ypQ dr J '' p"+y: 



(A70) 



where G approximates a cell-centered gradient from nodal data. Again, the f su- 
perscript refers to the fact that the updated velocity does not satisfy the divergence 
constraint. 

Finally, we use an approximate nodal projection to define U'^"''^ from U"'"'"^'^, such that 
U""*"^ approximately satisfies 

V • (/3o+'/^U"+i) = /3o+'/^ (^"+1 - :S^) , (A71) 

where /Sg is defined as 



C''' = ^°Y° ' ^^^' = ^ {Po^'^Po"-', rr\ g-^') . (A72) 

As part of the projection we also define the new-time perturbational pressure, Tr""'"/^ 
This projection necessarily differs from the MAC projection used in Step 3 and Step 
7 because the velocities in those steps are defined on edges and U"'"''^ is defined at cell 
centers, requiring different divergence and gradient operators. Details of the approxi- 
mate projection are given in Paper III. 

Step 12. Compute a new At. 

Compute At for the next time step with the procedure described in §3.4 of Paper III 
using wq as computed in Step 6 and U"^^ as computed in Step 11. 
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This completes one step of the algorithm. 

Figure [T^ is a flow chart summarizing the 12 step algorithm, including the notation used 
as we advance the solution by At. Figure [TS] is a flow chart of the advection steps (Steps 4 
and 8), which includes the notation we use as we advect the solution through a time interval 
of At. 



A. 5. Numerical Cutoffs 

As discussed in Paper IV, in order to prevent the velocity from becoming too large in 
low density regions far from the center of the star, we impose a cutoff at a moderately small 
density, Pcutoff, and hold the density at this constant value outside of the star. The cutoff 
affects the evolution in the following ways: 



After advancing enthalpy in the advection step (Steps 4Hii and 8Hii in §A.4p . we 
recompute the enthalpy using the equation of state if p < Pcutoff- 

When computing gravity ( §A.2p we only add po to rricnci if Po > Pcutoff in order to 
prevent an unphysical amount of mass from contributing to the calculation. 



When computing po in Enforce HSE ( §A.3.3p . we hold pq constant once po < PcutoS- 



In React State ( §A.3.ip . we set Uk = and pHnuc = if p < Pcutoff- 



When computing ip, (Steps 4F and 8F in §A.4p we set ip = if po < p, 



'cutoff • 



» When we compute the velocity forcing in Steps 3, 7, 11, and 12, we set the buoyancy 
term (the term proportional to p — po) to zero if p < Spcutoff- 

Additionally, we use an anelastic cutoff density, Paneiastic, in the computation of Pq (Steps 



3, 7, and 11 in ^A.4p. When poj < Panciastic, we set Po,j = {po,j/po,j~i)Po,j-i- 



B. Computing wq for spherical problems 



Recall that we want to solve 



1 "9 . 2 \ -n 1 f9po , dpo 



r^*(^^^'°)^^-ri;^lir^'"°i')- <^^' 



•</- 
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for the base state velocity, wq. We first decompose wq by setting wq = Wq + Swq, where the 
Wq term is the contribution to Wq due to the expansion term: 



19/ 2 \ T7in 

-IJTT [1" Wo) = S . 

iy-> ^ /^flY' ' 

Then we can write an equation for the remaining term, bw^. 



hl^ ('-^^-o) - -= 



TlPo [ dt 



dpo , .dpo' 



dr 



(B2) 



(B3) 



Multiplying equation fIBSP through by Tipo, taking another derivative with respect to r, and 
switching the order of temporal and spatial derivatives, we get: 



d_ 

dr 



^iPo 9,2. N 
^^ — (r 6wo) 



d dpo 
dt dr 



d_ 
dr 



[wq + 5wq 



dpo 
dr 



(B4) 



To solve for Swq we will need to substitute for the derivatives of po- To do so we start with 
the hydrostatic equilibrium equation. 



dpo 
dr 



-po9; 



9 



Gm, 



end 



(B5) 



where rrienciif') is the mass enclosed at radius r and G is the gravitational constant. Using 
this, we can then write equation ( ]B4p as: 



d 

dr 


TiPo d 2j: N 

[r^ Sr^^H 


d 
dt 

= 9 



{po9) + TT {woPo9) 



dr 



dpo 
dt 



d_ 
dr 



(wopo) 



Po 



dg , dg 



dt 



Wo 



dr 



(B6) 



The mass enclosed inside any radius, r, is rriendij') = Att J^ po{s)s'^ds, or alternately, dniend/dr 
Anr^po. The Lagrangian derivative of the enclosed mass is then: 

Domend 



Dt 



druend , drriend 



dt 

4vr — 

yotjo 



An 



dr 
Po{s)s^ds + Wor^Po 



A-K 



1 d{s'^poWo) 1 d{s'^rip) 



2 ds s^ ds 

= An {-s^poWo\l - s^?7p|[j + Wor^po) = -47rr\ 

where we used the spherical form of equation (29) in Paper III, 

dpo , 1 d{r^poWo) 1 d{r^r]p) _ 



—-s ds + wor po 
dt 



s^ds + wor^po 



dt 



+ 



dr 



dr 



(B7) 



(B8) 
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to eliminate dpo/dt. We note that in the absence of any mixing, rip = 0, and Donicnci/Dt = 0. 
Equation (IB7p allows us to write the Lagrangian change in the gravitational acceleration as: 



DQg__dg_ dg_ _ Dp f Gm^nd 

Dt ~ dt °9r ~ Dt\ r2 

2wpGmcnci 



Dq [ l\ G Doruenci 



- 4:TTGr]p 



2wog 



A-nGrjp. 



(B9) 



Putting it all together, equation (IB6P becomes 



d_ 

dr 






dpo . d 



Po 



-2wog 



A-nGr]^ 



Finally, we can use equation f lB8|) to write equation (IBlOp as: 
TiPo d 



d_ 
dr 



y,2 Qj, 



(r^^wo) 



9 






+Po 1 47rG?7p 



5' ^('^^^p) 4(wo + Swo)pog 



We discretize this elliptic equation in the radial dimension as: 



- 4:7rGpor]p. 



1 
Ar 



TiPod{r'^5wo)] \TiPod{r'^5wo)] 1 r4(r25ii'o)Po5' 



^i-y2 



r^ 1/ Ar - 

J -72 



dr 
[r'^Tlp),- (r\). 



c^r 



i-i. 



i-y2 



i-i 



(BIO) 



(Bii: 



(^) -(4.Gpo..),_,, (B12) 



where we choose to solve for [r'^dwo) rather than 5wq so that we can easily enforce d{r'^6wQ)/dr 
at the the upper boundary. Then, using hydrostatic equilibrium, we expand this to 

Ar 1 \ r^ / . Ar 






4 Po,j - Po,j-i 



S--y2 
^j-y2 

r^ 1, Ar 

j-y2 



Ar 



(r'^Sw, 



oJj-y2 



4 Poj - P0,j-1 \ ,2^ 



'j~% 



Ar 



r Wo)j^ii, 



(r^Vp)^ - {r^Vp)^_i - i^T^GpoVp)j-y, ■ 



(B13) 
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If we write this in matrix form, so that: 



Aj{r^6wo)j_:^/^ + Bj{r^6wo)j.i/^ + Cj{r^6wo)j+i/^ = Fj, 



then: 



A, 



B, 



C, 



F, 



1 /T^™ in 
1 / Ti p[f 



/\ ly Zi \ ly Zi 



j-i 



r~™ in 



r~"^ in 
Tl Po 



Ar2 

1 (rrp^ 



i-1 



4 Poj-Po,j^ i 

J -72 



Ar^ \ r^ 



4 Po!, -PoVi 
r^ 1, Ar 

. ^-72 



(r wojj-y. 



^j-y2 

r^ ,, Ar - 

J -72 



P<),-(-V),_i 



(B14) 

(B15) 
(B16) 
(B17) 



(B18) 



We define the lower boundary condition, 6wq = at r = 0, which corresponds to j = 0, by 
setting: 

Ao = Co = Fo = 0; Bo = 1. (B19) 

We also specify d{r'^6wo)/dr = at the the upper boundary, which corresponds to the 
location where po falls below Pcutos, by setting: 

An = -1; B^ = 1; Cn = Fn = 0. (B20) 

Finally, Wq^^ = Wo + Swq. Once po falls below Pcutoff, we hold r^ii^g"* constant. 



C. Test Problem Initial Model 

We use the same general initial model for the convergence test ( §6.3p . adaptive bubble 
rise test ( §6.4p and forced convection test ( §6.5p . We define a base temperature, Tbase = 
6 X 10® K, and density, pbase = 2.6 x 10^ g cm~^, at some height rbase above the bottom 
of the domain (rbase varies in each problem, and can be equal to zero). This defines a base 
entropy, Sbasc = ^(pbasc^basc)- The gravitational acceleration, g = —1.5 x 10^° cm s~^, is 
constant. The composition is uniform everywhere with X(^^C) = 0.3 and X(^^0) = 0.7. 
For r > rbase; we integrate the equation of hydrostatic equilibrium along with the constraint 
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that entropy is constant: 



1 

2' 



Po,j+i = Po,j + 7;^r{poj + poj+i)g. (CI) 



upwards from rbase- We define a temperature cutoff, Tcutoff = 10^ K, and when Tqj+i < T^utoff , 
we replace equation flC2p with Toj+i = Tcutos- If we choose rbase > 0, then we must also 
define the model for < r < rbase- In this case, we define a desired linear entropy profile 
with a discontinuous jump from Sbase as 

J- "^ '"base '5base /r^Q\ 

■^want TT'^base ~r ■ l^*^/ 

"J '"base J-^ 

This entropy profile creates a convectively stable layer below the atmosphere to prevent any 
motions generated from heating above from interfering with the lower boundary. The initial 
model for this region is then computed by integrating: 

Po,j = Po,j+i - 2^^(Po,j + Po,j+i)9, (C4) 

downward from r = rbase- 



-47- 



0.0x10 



0.0x10 



0.0x10 



0.0x10' 



1.0x10° 
radius (cm) 



I.Ox 10 
radius (cm) 



2.0x10° 




E 



2.0x10° 




a 
o 
a 

E 



2.0x10° 



2.6 xlO"^ 



2.5x10^ 



2.4x10' 



t=0 

t=2, 1 D MAESTRO 

t=2, 1 D CASTRO 

t=2, 3D MAESTRO 



0.0x10 




7.5x10° 
radius (cm) 



1.5x10 



1.7x10 



1.6x10" 



t=0 ■ 

t=2, 1 D MAESTRO 

t=2, 1 D CASTRO □ 
t=2, 3D MAESTRO x 




0.0x10 



7.5x10° 
radius (cm) 



1.5x10 



9.0x10° 



8.0x10° 



7.0 X 10° 



6.0x10° 



5.0x10° 



t=0 ■ 

t=2, 1 D MAESTRO 

t=2, 1 D CASTRO □ 
: t=2, 3D MAESTRO x 




3.0x10' 



radius (cm) 



Fig. 7. — Plots of po (Top), pq (Middle), and T (Bottom) vs. radius for a white dwarf star 
subject to heating. The initial profiles are on the left. Close-up views of the initial profiles 
and final solutions are on the right. We use three test codes: a one-dimensional version 
of MAESTRO in spherical coordinates, a one-dimensional version of the compressible code, 
CASTRO, in spherical coordinates, and a full-star three-dimensional version of MAESTRO. 



-48- 





Fig. 8. — Profile of T — T for a fiot bubble in a white dwarf environment. The black and 
red lines represent grids of increasing refinement. The vertical distance shown is from z = 
to 9.2 X 10^ cm. 




^M 








q 


^ 






^! 


V 













r 


Hu.^^.J 








ta. J 



Fig. 9. — Time-lapse cross section of Figure [H] at t = 0,0.5,1.0,1.5,2.0, and 2.5 s for a 
single-level simulation (above) an adaptive simulation with two levels of refinement at the 
same effective resolution (below). 
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Fig. 10. — Comparison between po (Top), po (Middle), and T (Bottom) at t = 0, and t = 30 s 
in a white dwarf environment with a heating layer for the single-level fine grid and multi-level 
simulations. The multi-level simulation captures the same expansion of the base state as the 
single-level simulation. 
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Fig. 11. — Temperature plots at t = 3 s (above) and t = 4 s (below) of a white dwarf 
environment with a heating layer. The single-level coarse grid solution is on the left, the 
multi- level solution is in the center, and the single- level fine grid solution as on the right. The 
colored boxes indicate the grid structure at each level of refinement. The vertical distance 
shown is from z = 5 x 10^ cm to 2.2 x 10^ cm. At t = 3 s, the multi-level simulation is 
able to capture the finer-scale structure visible in the single-level fine grid simulation, which 
is not captured as well in the single-level coarse grid simulation. At t = 4 s, a finer-scale 
structure is still visible in the multi-level simulation, but the solution begins to drift from 
the single-level fine grid simulation, which is expected since we are deliberately refining only 
a part of the convective region. 
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Fig. 12. — Initial grid structure for a full white dwarf star simulation with 2 levels of 
refinement. The colored boxes indicated the grid structure at each level of refinement, each 
grid containing up to 32^ cells. The finest grids have effective resolution of 384^. The shaded 
sphere indicates where the density is 10^ g cm~^ or greater. In this simulation, we have 
chosen to include all points where p > 10^ g cm~^ at the first level of refinement and all 
points were where p > 10^ g cm~'^ at the second level of refinement. There are 216 black 
grids at the coarse level, 125 red grids at the next level, and 212 blue grids at the finest level. 
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Fig. 13. — Peak temperature (Top) and peak radial velocity (Bottom) as a function of time 
for a single-level and adaptive simulation, each with an effective 384^ resolution. 
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Fig. 14. — A flowchart of tlie algoritliiii. Tlie tliermodynainic state variables, base state 
variables, and local velocity are indicated in each step. Red text indicates that quantity was 
updated during that step. The predictor-corrector steps are outlined by the dotted box. The 
blue text indicates state variables that are the same in Step 6 as they are in Step 2, i.e., 
they are unchanged by the predictor steps. 



-54- 



StartofStep4/8 

p">, (phf>, (pXf, V" 



I 



A. Advect Base Density 

p">. {phf>, (pXf>. V" 



I 



B. Update Species 

p"', iphf, (pXf, t" 



I 



C. Define q 

p<'>, (phf, (pxf>, r" 



D. Correct Base Density 

p'", (phf>, (pXf>, V" 

P:"' Pn"' (Pf'tn" 



E. Enforce HSE 

p">, (phf, (pxf>, r" 

Po"' Po"' (P^h" 



I 



F. Compute ijj 

p'^ (ph)">, (pxf>, r" 



, n+1 „ n+1 



I 



G. Advect Base Enthalpy 

p», (phf, (pxf>, r" 

Pn""- Pn""' (Ph)n"' 



I 



H. Update Enthalpy 

p'", (phf, (pxf>, v" 

P:"- Pn""' (pK" 



i 



I. Update T 

p'", (phf, (pXf>, T<'> 

p:", Pn"". (Ph):" 



Fig. 15. — A flowchart for Steps 4 and 8. The thermodjTiainic state variables and base 
state variables are indicated in each step. Red text indicates that quantity was updated 
during that step. Note, for Step 4, the updated quantities should also have a -k superscript, 
e.g.. Step 81 defines T^^) while Step 41 defines T^^)'* . 
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